R Notebook

library(Seurat)

mat <- list()
mat$HC1 <- Read10X_h5(filename = "../data/COVID-19/GSM4475048_C51_filtered_feature_bc_matrix.h5")
mat$HC2 <- Read10X_h5(filename = "../data/COVID-19/GSM4475049_C52_filtered_feature_bc_matrix.h5")
mat$HC3 <- Read10X_h5(filename = "../data/COVID-19/GSM4475050_C100_filtered_feature_bc_matrix.h5")
mat$HC4 <- Read10X(data.dir = "../data/COVID-19/GSM3660650/")

mat$M1 <- Read10X_h5(filename = "../data/COVID-19/GSM4339769_C141_filtered_feature_bc_matrix.h5")
mat$M2 <- Read10X_h5(filename = "../data/COVID-19/GSM4339770_C142_filtered_feature_bc_matrix.h5")
mat$M3 <- Read10X_h5(filename = "../data/COVID-19/GSM4339772_C144_filtered_feature_bc_matrix.h5")

mat$S1 <- Read10X_h5(filename = "../data/COVID-19/GSM4339773_C145_filtered_feature_bc_matrix.h5")
mat$S2 <- Read10X_h5(filename = "../data/COVID-19/GSM4339771_C143_filtered_feature_bc_matrix.h5")
mat$S3 <- Read10X_h5(filename = "../data/COVID-19/GSM4339774_C146_filtered_feature_bc_matrix.h5")
mat$S4 <- Read10X_h5(filename = "../data/COVID-19/GSM4475051_C148_filtered_feature_bc_matrix.h5")
mat$S5 <- Read10X_h5(filename = "../data/COVID-19/GSM4475052_C149_filtered_feature_bc_matrix.h5")
mat$S6 <- Read10X_h5(filename = "../data/COVID-19/GSM4475053_C152_filtered_feature_bc_matrix.h5")

meta.data <- read.delim("../data/COVID-19/all.cell.annotation.meta.txt")
rownames(meta.data) <- meta.data$ID
head(meta.data)

numbering = c(HC1 = 1,
              HC2 = 2,
              HC3 = 3,
              HC4 = 4,
              M1 = 5,
              M2 = 6,
              M3 = 7,
              S1 = 9,
              S2 = 8,
              S3 = 10,
              S4 = 11,
              S5 = 12,
              S6 = 13)

for (i in names(mat)){
  colnames(mat[[i]]) <- gsub('-', '_', colnames(mat[[i]]))
  colnames(mat[[i]]) <- gsub('1', numbering[i], colnames(mat[[i]]))
}

for (i in names(mat)){
 mat[[i]] <- mat[[i]][, meta.data$ID[meta.data$sample_new == i]]
}

objs <- list()
for (i in names(mat)){
 objs[[i]] <- CreateSeuratObject(mat[[i]], project = i, meta.data = meta.data[meta.data$sample_new == i, ])
}

obj <- merge(objs[[1]], objs[-1])

obj <- NormalizeData(obj, verbose=FALSE)
obj <- FindVariableFeatures(obj, selection.method = "vst", nfeatures = 2000, verbose=FALSE)

obj <- ScaleData(obj, verbose=FALSE)
obj <- RunPCA(obj, features = VariableFeatures(object = obj), verbose=FALSE)

ptm = proc.time()
obj <- RunHarmony(obj, "sample_new", plot_convergence = TRUE)

Harmony 1/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 2/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 3/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony converged after 3 iterations

proc.time() - ptm

   user  system elapsed 
129.656  12.969 143.064 

obj <- RunUMAP(obj, reduction = "harmony", dims = 1:30, verbose=FALSE, reduction.name = "umap", reduction.key = "ALTUMAP_")

Cannot add objects with duplicate keys (offending key: ALTUMAP_), setting key to 'umap_'

DimPlot(obj, reduction = "umap", group.by = "sample_new", label.size = 5,  pt.size = 0) +
  theme(axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank(),
        text=element_text(size=18))
ggsave("covid_harmony_sample.pdf")

Saving 7.29 x 4.5 in image

DimPlot(obj, reduction = "umap", group.by = "sample_new", label.size = 10,  pt.size = 0) +
  theme(axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank(),
        text=element_text(size=30))
ggsave("covid_harmony_sample.png", height=12, width=18, dpi = 100)

DimPlot(obj, reduction = "umap", group.by = "celltype", label = T, repel=T, label.size = 5,  pt.size = 0) +
  theme(axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank(),
        text=element_text(size=18))

Using `as.character()` on a quosure is deprecated as of rlang 0.3.0.
Please use `as_label()` or `as_name()` instead.
This warning is displayed once per session.

ggsave("covid_harmony_label.pdf")

Saving 7.29 x 4.5 in image

DimPlot(obj, reduction = "umap", group.by = "celltype", label = T, repel=T, label.size = 10,  pt.size = 0) +
  theme(axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank(),
        text=element_text(size=30))
ggsave("covid_harmony_label.png", height=12, width=18, dpi = 100)